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Response of a circumbinary accretion disc to black hole mass loss 
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ABSTRACT 

We investigate the evolution of the surface density of a circumbinary accretion disc after the 
mass loss induced by the merger of two supermassive black holes. We first introduce an ana- 
lytical model, under the assumption of a disc composed of test particles, to derive the surface 
density evolution of the disc following the mass loss. The model predicts the formation of 
sharp density peaks in the disc; the model also allows us to compute the typical timescale for 
the formation of these peaks. To test and validate the model, we run numerical simulations of 
the process using the Smoothed Particle Hydrodynamics (SPH) code PHANTOM, taking fluid 
effects into account. We find good agreement in the shape and position of the peaks between 
the model and the simulations. In a fluid disc, however, the epicyclic oscillations induced by 
the mass loss can dissipate, and only some of the predicted peaks form in the simulation. 
To quantify how fast this dissipation proceeds, we introduce an appropriate parameter, and 
we show that it is effective in explaining the differences between the analytical, collisionless 
model and a real fluid disc. 

Key words: accretion, accretion discs - black hole physics - hydrodynamics. 



1 INTRODUCTION 

Supermassive black holes (SMBHs) are hosted in the nuclei of 
most galaxies (Kormendy & Richstone 1995). If two host galax- 
ies merge, as predicted by hierarchical galaxy formation models, 
the two black holes can form a supermassive black hole binary 
(Begelman et al. 1980) via dynamical friction. Mergers of gas rich 
galaxies also drive a large quantity of gas in the centre, potentially 
forming massive, rotationally supported disc (Escala et al. 2005; 
Dotti et al. 2007). If this gas can cool efficiently, it settles into a 
thin accretion disc surrounding the black hole binary. 

If the binary hardens up to sub-parsec scales through gravi- 
tational slingshot ejection of stars (Milosavljevic & Merritt 2003) 
or other processes, such as the interactions with a gaseous accre- 
tion disc (Ivanov et al. 1999; Lodato et al. 2009; Nixon et al. 201 1), 
eventually gravitational waves (GW; Peters 1964) will carry away 
the remaining orbital energy and induce the coalescence in less than 
a Hubble time. The detection of these GWs is expected by pro- 
posed experiments such as the Laser Interferometer Space Antenna 
(LISA). 

General relativity predicts (Pretorius 2005; Campanelli et al. 
2006; Tichy & Marronetti 2008) that the gravitational waves emit- 
ted during the merger of the black holes carry away energy and 
momentum, thus resulting in a mass loss and in a recoil of the 
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remnant black hole. Many authors (e.g. Schnittman & Krolik 2008; 
Rossi et al. 2010; Corrales et al. 2010; Zanotti et al. 2010) have ex- 
plored the second possibility, with the final aim of predicting an 
electromagnetic (EM) afterglow of the coalescence. The detection 
of such an afterglow would help in constraining the properties of 
the merged black holes and of the host galaxies. In this paper 
we explore the first possibility, studying the response of the ac- 
cretion disc to the mass loss. For the rest of this paper, we ne- 
glect the effect of recoil. Previous work for this physical case 
have been done by Megevand et al. (2009); O'Neill et al. (2009); 
Corrales et al. (2010) by means of numerical simulations. In the 
present paper, we present an analytical model for the surface den- 
sity evolution of the disc following the merger, derived under the 
assumption of a collisionless disc. To assess the validity of our 
model, we compare it with the outcome of 3D numerical hydrody- 
namical simulations, using the Smoothed Particle Hydrodynamics 
(SPH) code PHANTOM (Price & Federrath 2010; Lodato & Price 
2010). 

The geometry of the accretion disc before the merger 
has been discussed by Armitage & Natarajan (2002) and by 
Milosavljevic & Phinney (2005). We assume that the plane 
of the disc coincides with the orbital plane of the binary 
(Larwood & Papaloizou 1997; Ivanov et al. 1999). For such a 
configuration, the secondary black hole will open a gap 
(Artymowicz & Lubow 1994) in the disc, creating a hollowed re- 
gion of size approximately 2a, where a is the semi-major axis of 
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the binary. The evolution of the system will then progress at the 
viscous timescale, the disc and the binary being in contact due to 
tidal interactions. When the gravitational wave emission sets in, the 
structure of the circumbinary disc and of the binary rapidly decou- 
ple, and the coalescence occurs on a very fast timescale compared 
to the disc dynamical one. We can then assume that at the moment 
immediately preceding the merger the disc is in rotation around a 
point mass, whose magnitude is given by the total mass of the bi- 
nary. For our purposes, the merger of the two black holes can be 
described as an instantaneous reduction of the mass of the central 
object. 

The paper is organised as follows. In Section 2, we derive an 
analytical model for the evolution of the surface density of the disc 
after the merger, and thus the mass loss, occurs. The model is then 
compared with numerical simulations, described in Section 3; we 
present our results in section 4. Finally, in Section 5, we draw our 
conclusions. 



2 ANALYTICAL MODEL 

In this section we develop a ID analytical model for the surface 
density evolution of the accretion disc in the case of mass loss. We 
make the assumption that the disc is composed of test particles, 
and we neglect the effects of pressure; thus we deal with a colli- 
sionless disc. This is the same approach that Lippai et al. (2008) 
followed for analysing the response to a black hole recoil (in that 
case through the use of numerical simulations). Since real discs are 
collisional, the validity of our model has to be investigated and will 
be discussed in section 4. 



2.1 Derivation 

Before the mass loss happens, we make the simplifying hypothesis 
that the test particles in the disc are on circular orbits. That is, we 
assume a Keplerian rotation curve with a rotation speed given by 
V K = y/GM/R, where G is the gravitational constant, M the 
total mass of the two merging black holes and R the distance of 
the particle from the center of mass of the black holes in the disc 
plane. The specific angular momentum j of each particle is given 
VkR = V GMR, or, conversely, the orbital radius as a 



by j 

function of the specific angular momentum reads 

•2 



R = 



GM' 



(1) 



The merger of the two black holes, and thus the mass loss, 
happens at t = 0, and causes a fractional change in the mass of the 
black holes of magnitude e <C 1. Thus for t > the new mass of 
the remnant M' is given by M' = M(l — e). 

The mass loss does not change the specific angular momen- 
tum of the particles, as happens in the case where a recoil is present 
(Rossi et al. 2010). However, due to the change in mass, every par- 
ticle initially at radius R' will have a circularisation radius ex- 
pressed by: 



RciTc{R') 



GM(1 



R'(l + e) 



(2) 



Thus, coherently with a decrease in the central mass, the cir- 
cularisation radius increases by the same fractional amount, so that 
after the merger the particles are in elliptical orbits. To first order in 
e, these orbits can be described using the epicycle approximation. 
This means that the radius oscillates in an harmonic fashion about 



the new circularisation radius with the epicyclic frequency k. The 
distance of the particle from the center can be then expressed as a 
function of the initial radius R' and of time as: 



Rnew(R' ,t) — Rcirc(R') + 

A(R cilc (R'))sm(K.t + 3tt/2), 



(3) 



where A(Rciic(R )) is the amplitude of the epicyclic oscillations. 
The initial phase is required to ensure that at t = the particle is at 
R = R' and its radial velocity vanishes. These requirements also 
sets the amplitude of the oscillations, that is: 



A — R c i IC — R' — tRcir 



(4) 



Finally, we recall that for a Keplerian disc the epicyclic fre- 
quency is equal to the orbital frequency Q. Evaluating this quantity 
at the circularisation radius one obtains: 



GM(1 



(5) 



We have now all the ingredients to compute how the mass loss 
changes the surface density E of the disc given the initial surface 
density Eo(-R). Under the constraint that the mass of test particles 
does not change, we can evaluate the required quantity as follows: 



T,(R,t) = 



flout 



Xo(R')S[R- R ncv {R' ,t)]dR' , 



(6) 



where the limits of the integral are the inner iZi n and the outer -R ou t 
radius of the disc, S is the Dirac function and the extra R' /R factor 
accounts for the surface element variation and enforces the conser- 
vation of mass. The Dirac function in the integral takes care of se- 
lecting the contributions to the surface density of particles, initially 
at radius R', that at time t are at radius R. 



2.2 Implied timescale 

There is no analytic solution for the inversion of the argument of 
the Dirac delta function. However, we show that our model predicts 
the formation of density peaks in the disc. This is a consequence of 
the strong radial dependence of the epicyclic frequency with ra- 
dius, k oc i?~ 3 / 2 . This implies that, given two particles at a certain 
distance AR, they will quickly go out of phase. When the phase 
difference between two given particles is of order n, one particle 
will be increasing its radius while the other one is diminishing it. If 
they are close enough, they will intersect, forming a region of high 
density. Conservation of mass implies that we expect also regions 
of depression. 

With a simple argument we can derive a timescale t^p for the 
formation of these density peaks. The formation of a peak occurs 
when the phase difference A<j> is of the order n for particles that 
are 2ei? c i rc ~ 2eR apart, that is, the maximum initial distance for 
which they can intersect. Then we can write: 



A<j} = t dp [n(R) - k(R + 2eR)] = 



dyn 



(7) 



where ta yn = I/O is the dynamical timescale. Solving for td P we 
obtain, neglecting factors of order unity: 

tdyn 



tdp 



(8) 



This is roughly the same time scale found by Schnittman & Krolik 
(2008); however we stress that in this section we built a time- 
dependent model that gives the time evolution of the density and 
not only the timescale for the formation of density peaks. 
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used 1CP J in our computations). Once a has been chosen, the inte- 
gral can be computed with standard numerical techniques (we used 
Simpson's rule). To ensure the discretisation does not significantly 
modify the final result, we used a sufficiently large number of sam- 
pling points such that the separation between them is much less 
than a (we used a 1CP 2 ratio in our computations). 

Figure 1 shows the numerical solution, normalized to the ini- 
tial unperturbed value, at time t = 48d after the merger (red solid 
line), for fiducial parameters of M = 10 6 A/q and e = 0.05. The 
unperturbed surface density profile is taken to be E(_R) oc R~ p . 
As shown in 2.2, sharp over-density peaks develop in the disc, as 
well as region of depression. The figure also shows a resolution 
study changing the r parameter that controls the Gaussian width a: 
dashed line green is the result for r = 10~ 4 and blue dotted line for 
t = 10~ 3 . Outside 200 R s we find excellent convergence, while 
some difference remains at the smallest radii. 



Figure 1. Numerical solution, normalized to the initial unperturbed value, 
at time t = A8d after the merger (red solid line). Sharp over-density peaks 
develop in the disc, as well as region of depression. The figure also shows a 
resolution study changing the r parameter that controls the Gaussian width 
a: dashed line green is the result for r = 10~ 4 and blue dotted line 
for t = 10 -3 , while the red line is the result for our fiducial value for 
r = 1CT 5 . Outside 200 R s we find excellent convergence, while some 
difference remains at the smallest radii. 



To get an estimate for this timescale, we can use the 
Armitage & Natarajan (2002) model for the decoupling separation 
of the binary a, that corresponds roughly to the inner radius of the 
disc. Their expression is however valid only for the case of extreme 
mass ratio binary. Using their arguments, one can derive the fol- 
lowing expression, that has a general validity: 



(—) 

\15V2J 



2/5 



-2/5 



-4/5 
17 J 1 



2/5 



(1+9) 



-4/52GM 



Here a is the Shakura & Sunyaev (1973) parameter for viscosity, 
H/R is the aspect ratio of the disc prior to decoupling and q is the 
mass ratio of the binary: q — M2/M1, where Mi and M2 are the 
masses of the two black holes. 

For fiducial parameters of a = 10~ 2 , H/R = 10 -2 , M = 
10 6 M Q and q = 1/2, we obtain an inner disc radius of ~ 2 AU, 
so that we expect td P — 3 days at the inner disc radius for 
e = 0.05. Note that in units of 2GM/c 2 , that is approximately 
the Schwarzschild radius R s of the remnant, a ~ 100i? s for a large 
range of values of q, due to its weak dependence on q. 



2.3 Numerical solution 

Since the integral (6) has no analytic solutions, we computed it 
numerically. We describe here briefly the method we used. To deal 
with the Dirac delta function, we replaced it using a Gaussian, so 
that the integral becomes: 



T,(R,t) = 



1 



exp 



27TCT 

(R- 



-Rill 

Rncv 



2a 2 



dR', 



(10) 



where a is the Gaussian width. Since the actual integral is recov- 
ered in the limit in which a — > 0, a is chosen to be much smaller 
than the radius R at which we are interested in computing the sur- 
face density: rr = tR, where r is an arbitrary small number (we 



3 SIMULATION DETAILS 

In order to test the validity of the model developed in the previous 
section, we set up numerical hydrodynamic simulations. Simula- 
tions are capable of following the full 3D, possibly non-linear hy- 
drodynamical evolution of the disc, and can then be used to quan- 
tify the limits of a collisionless model. 

3.1 PHANTOM 

Our simulation use the PHANTOM SPH code. SPH 
(Gingold & Monaghan 1977; Benz 1990; Monaghan 2005; 
Price 2012) is a well-known Lagrangian technique to solve 
the hydrodynamics equations. It is based on a discretisation of 
the fluid mass in terms of particles, where each particle rep- 
resents a smeared out distribution of density. The formulation 
of SPH implemented in PHANTOM is the so-called 'grad-h' 
formulation (Price & Monaghan 2004, 2007), which is based on 
a self-consistent derivation of the equations of motion from a 
least action principle. This guarantees an exact conservation of 
momentum, angular momentum, and energy. In this formulation 
the smoothing length h, which sets the effective resolution length 
of the simulation, and the density p are mutually dependent, so 
that the smoothing length is adaptively adjusted to give a better 
resolution in the high density regions. 

PHANTOM (Price & Federrath 2010; Lodato & Price 2010) is 
a low-memory, highly efficient 3D SPH code specifically designed 
for non self-gravitating problems.The code implements a linked list 
(Monaghan & Lattanzio 1985) for neighbour finding. The code has 
been parallelised using both OpenMP and MPI, although in this 
paper only the shared memory OpenMP parallelisation has been 
used. 



3.2 Initial conditions 

Our initial conditions consist of a thin disc (Pringle 1981). Due to 
the arguments discussed in section 2.2, we set the inner radius of 
the disc to 100i? s , where R s is the Schwarzschild radius of the 
black hole. We set the outer radius to be 10 4 i? s . 

The radial distributions of surface density and sound speed, 
c s , are power laws: 



E(i?) oc R~ p , 
Cs(R) oc R~", 



(11) 
(12) 
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Figure 2. Face-on view of the disc for simulation cold-iso at four different times. The rendered quantity is surface density, normalized to a total disc mass of 
500 Mq . The formation of sharp density peaks is particularly evident. 



Abbreviation 


flout 


No. of particles 


H/R 


EOS 


R/h 


H/h 




cold-iso 


10 4 


2M 


0.032 


Isothermal 


~ 60 


2 


0.025 


cold-iso-hres 


10 3 


2M 


0.032 


Isothermal 


~ 300 


10 


0.005 


hot-iso 


10 4 


7M 


0.125 


Isothermal 


~ 60 


8 


0.00625 


hot-iso-hres 


10 3 


7M 


0.125 


Isothermal 


~ 300 


37 


0.00125 


cold-adiabatic 


10 4 


2M 


0.032 


Adiabatic 


~ 60 


2 


0.025 


hot-adiabatic 


10 4 


7M 


0.125 


Adiabatic 


~ 60 


g 


0.00625 



Table 1. The table summarizes the simulations run and the abbreviation used in the text to refer to them. We explored different temperatures of the disc and 
different equations of state. To check for numerical convergence, we increased radial resolution of our simulations by truncating the disc at the radius -Rout 
indicated (in units of the Schwarzschild radius). The aspect ratio H/R and the radial resolution R/h are evaluated indicated in the table are evaluated at 10 2 i?, s , 
while the vertical resolution H/h is indipendent of radius. Note that the radial resolution is an increasing function of radius, so that the value reported in the 
table is a lower limit for this value. For details about the equation of state, see section 3.2; about the radial resolution, section 4.2. We report in the table also 
the values of dgg Bax , the maximum equivalent Shakura-Sunyaev parameter that is attained at shocks. The typical values of the effective ctgg away from 
shocks are one or two orders of magnitude smaller. 



where the exponents p and q are chosen to satisfy the relation p + 
2q — 3. This ensures that the vertical averaged smoothing length 
(h) oc H, where H is the height of the disc. This ensures that 
the disc is equally vertically resolved at each radius. For all our 
simulations we assume p — 3/2 and therefore q = 3/4. We will 
discuss in section 4 the normalisation of the sound speed chosen. 
The disc is isothermal in the vertical direction. Initial conditions 
were generated using random placement of particles. 

We assume a Keplerian rotation curve, taking into ac- 
count also the pressure gradient corrections, and we neglect rel- 
ativistic corrections. A simple estimate using the well-known 
Paczyhsky & Wiita (1980) potential gives that the error is of the 
order of a percent for the innermost radius of the disc, and even 



smaller for the rest of the disc. Thus the use of a Newton potential 
for modelling the black hole is sufficient. 

For the thermodynamics, we experimented both with the use 
of an isothermal equation of state, where every particle keeps fixed 
its initial temperature given by equation (12), and the use of an adi- 
abatic equation of state, where the internal energy of each particle 
can change. In both cases we assumed an ideal gas law. The two 
equations of state chosen correspond to the physical case of a very 
fast cooling timescale and a very long one. We choose the value 
7 = 5/3, corresponding to ideal monatomic gas. 

At time t = 0, we assume that the central source of the gravi- 
tational field experiences a mass loss of fractional change e = 0.05, 
and we let the system evolve in time. In section 2.2 we have shown 
that the timescales we are interested in are of the order of days, 
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Figure 3. Face-on view of the disc for simulation hot-iso at four different times. The rendered quantity is surface density, normalized to a total disc mass of 
500 Mq . 



thus much shorter than the viscous timescales. For this reason we 
do not include physical viscosity in our simulations, and we in- 
clude only the artificial viscosity needed to capture shocks in SPH. 
PHANTOM includes a slightly modified version of the Monaghan 
(1997) formulation of artificial viscosity, based on an analogy with 
Riemann solvers. Details on the implementation can be found in 
Price & Federrath (2010). We chose standard values to parametrize 
the viscosity, with a = 1 and j3 — 2. The meaning of a can 
be understood by comparison with the compressible Navier-Stokes 
equation. It can be shown that in the case of constant a the artifi- 
cial viscosity (e.g. Lodato & Price 2010) is equivalent to a Shakura- 
Sunyaev parameter ass, given by: 



ass = ^a-. 



(13) 



The corresponding bulk viscosity coefficient is a factor 5/3 greater 
than the shear. 

In our simulation we use the Morris & Monaghan (1997) 
switch to reduce dissipation away from shocks, so that the indi- 
vidual a of each particle can vary between 0.01 and 1, so that the 
viscosity coefficients vary. In a presence of a shock front, however, 
the values of a are near the maximum, so that one can use the max- 
imum value to compute the viscosity coefficients. 

We report in table 1 the values of ass, max for the various 
simulations, where ass, max is the maximum equivalent Shakura- 
Sunyaev parameter that is attained at shocks, where a = 1. Since 
in our simulations we do not have strong shocks, the typical values 
of the effective ass are one or two orders of magnitude smaller. 



4 RESULTS 

In this section we present our results and we discuss the compar- 
ison between our model and the simulations. We run simulations 
with four different sets of physical parameters, studying the depen- 
dence on the normalisation of the sound speed and on the equation 
of state. For the sound speed, we took two different values, that 
correspond roughly to an aspect ratio H/R of the disc of 0.032 
and 0.125 at R — 10 2 R B . In the rest of the disc, the aspect ratio 
exhibits a mild dependence on radius (oc R^ 1 ^ 4 ), given the chosen 
sound speed profile. We will refer in the rest of the paper to the two 
cases as "cold" and "hot" disc. 

We studied also the resolution dependence in our simulations 
to check that numerical convergence has been reached. To run sim- 
ulations at a significantly higher resolution than our "fiducial" one, 
we simulate only the inner part of the disc, truncating it at 10 3 i? s , 
keeping the same number of particles. This ensures that the mass 
of each SPH particle becomes significantly lower. We summarize 
the simulations run in table 1 . 

Since we employed Schwarzschild units for the lengths, our 
results do not depend on the mass of the central black hole. To show 
physically meaningful times, however, in the following section time 
and temperature are expressed in physical units for the fiducial case 
in which M = 10 6 M o . 

4.1 Isothermal simulations 

We first present results from the isothermal simulations, namely 
cold-iso and hot-iso. Figure 2 and Figure 3 present face-on views of 
the disc, rendering surface density, at four different times for these 



© 201 1 RAS, MNRAS 000, 1-10 




two simulations. More quantitatively. Figure 4 shows a compari- 
son between the two simulations (cold-iso: solid red line; hot-iso: 
dashed green line) and the analytic model (dotted blue line) at four 
different times, namely t=24d, 48d, 72 and 97 d after the mass loss. 
We show surface density, normalized to the pre-massloss value, as 
a function of radius. 

The main effect visible is the formation of density peaks after 
the mass loss. Once formed, these peaks then propagate outwards. 
This is consistent with the prediction of the model that the peaks 
first form in the inner part of the disc, where the timescale for their 
formation is faster (see equation 8), and then travel outwards. There 
is a quite remarkable agreement between the model and the simu- 
lations in the outer part of the disc. We should stress that, since 
our analytic model is collisionless, the propagation of these peaks 
is not a hydrodynamical effect (like the propagation of waves in a 
fluid medium). Instead, this effect comes only from the dynamics 
of the particles comprising the disc. However, fluid effects in a real 
disc can modify the evolution in time, and this is the reason why we 
need also hydrodynamical simulations in addition to the analytical 
model. Before proceeding with a more detailed analysis, we note 
that also Corrales et al. (2010) run hydrodynamics simulations in 
the case of mass-loss, with a very different numerical scheme (they 
use the grid-code FLASH). Qualitatively, the outcome of their sim- 
ulations (see their Figure 5) is very similar to what we have ob- 
tained. This a confirmation that the result does not depend on the 
numerical method employed. 

We find that a fluid disc differs from the model in two main 
aspects: 



(i) the position of the peaks for the hot disc does not coincide 
with the one predicted by the model (see figure 4); 

(ii) the disturbances are damped in the fluid disc, so that in the 
inner part of the disc the analytical model and the simulations dis- 
agree. This effect is stronger for the hot disc. 

To account for the first difference, we can notice, using equa- 
tion (8), that disturbances travel at a speed Vd given by: 

v d = feV K , (14) 

where / is a number of order unity. On the other hand, no infor- 
mation in a fluid can travel faster than the sound speed c s . It is 
then relevant to consider the ratio M e — eVk/c s . If this ratio is 
greater than 1, the disturbance travels faster than the sound speed. 
In this regime, the inner fluid cannot inform the outer one of the 
disturbance, and the position of the peaks is the one given by the 
analytical model. In the opposite regime, the fluid will communi- 
cate to the outer disc the formation of the peak, thus distorting its 
shape and shifting it to outer radii. 

In our discs, M c oc R' 1,2 /R~ 3/A oc R 1/4 , with a value 
at 10 3 i? s given by 2.8 for the cold case and 0.7 for the hot case. 
Figure 5 shows M t as function of radius for the two simulations. 
It can be readily seen that for the cold disc A4 t is always greater 
than 1, and this is consistent with what we have pictured above. For 
the hot disc, the results are intermediate between the two extremes, 
and this accounts for the shift of the peaks and for the distortion of 
the shape. 

In particular, if the velocity of the peak is given by equation 
(14), the position of the peak as a function of time is given by 
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Figure 5. The M ( parameter as a function of radius for the cold-iso (red 
solid line) and hot-iso (green dotted line) simulations. While for the cold 
disc the parameter is everywhere greater than 1, thus being in the regime 
in which the disturbance travels at high supersonic speed, in the hot case 
the parameter is slightly less than 1 , and the disturbance travels at a speed 
comparable with the speed of sound. 
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Figure 6. Position of the outermost peak as a function of time for the cold- 
iso simulation (crosses) and for hot-iso (circles). Overplotted are the best fit 
to the points using equation (15). In the cold case we get / = 1.66, while 
in the hot case we get a higher value of / = 1.87. The value of Ro we 
get from the fit is compatible in both cases, given the fit uncertainties, with 
zero. 



Rp(t) 



2/3 



(15) 



Figure 6 shows the position of the outermost peak as a function of 
time for the two simulations cold-iso (crosses) and hot-iso (circles). 
Overplotted are the best fit to the points using equation (15), where 
we let / and Rq be free parameters of the fit. In both cases there is 
a good agreement between the fit and the position of the peak. The 
value of / we get from the fit are 1.66 for cold-iso and 1.87 for hot- 
iso, that is, in the hot disc case the peaks moves at a slightly higher 
velocity. These values also confirm that the constant / is of order 
unity, and thus that td P , from which equation (14) has been derived, 
is a valid timescale for the formation of these density peaks. The 
value of Rq we get from the fit is compatible in both cases, given 
the fit uncertainties, with zero. 

Concerning the second difference, this is due to the fact that 
in the fluid disc the gas compression that occurs in the density 




Figure 7. Surface density normalized to initial value at time t = 24d for 
simulations cold-iso (solid red line), cold-iso-hres (dashed green line), hot- 
iso (blue dotted line) and hot-iso-hres (light blue dotted-dashed line). No 
significant difference is found increasing the resolution. 



peaks causes the energy of the epicyclic oscillations to be dissi- 
pated. Thus, for a given radius only a certain number of peaks can 
form before the oscillations are damped. This is due to a combina- 
tion of shock dissipation and to p dV work. While in the isothermal 
case this work is lost from the system, in the adiabatic case it is con- 
verted into heat. As the gas temperature increase, the pressure term 
becomes more and more important, so that in the adiabatic simula- 
tions the oscillations are damped much faster than in the isothermal 
case. The pressure term also explains why in the hot simulations the 
perturbations are damped faster than in the cold one. Since the dy- 
namical timescales are much faster in the inner disc, in the inner 
region the epicyclic motions already dissipated their energy, while 
in the outer region we can still see their effect. 



4.2 Resolution dependence 

A very important issue is establishing if our simulations have 
reached numerical convergence. Too low resolution would imply 
that we cannot resolve properly some of the density peaks, thus be- 
ing unable to distinguish between physical and numerical effects. 
We now turn to this problem. 

A helpful parameter to consider is the radial resolution R/h, 
where h is the SPH smoothing length. This is connected to the ver- 
tical resolution through the H/R aspect ratio, so that the radial 
resolution can be expressed as R/h = (H/R)~ H/h. It is then 
not possible, for a given aspect ratio, to have arbitrary values of 
the radial and vertical resolution. In this paper, when comparing 
simulations with different aspect ratios, we used the same radial 
resolution, so that the numerical uncertainties are the same for the 
two cases. Thus, the differences seen in the previous section be- 
tween the hot and cold cases cannot be due to a resolution problem. 
Table 1 also reports the radial resolutions of our simulations. No- 
tice that since the aspect ratio is a function of radius, the radial 
resolution depends on radius (while our setup ensures that H/h = 
const). However, the scaling is very mild and does not affect this 
discussion. In addition, we checked that all of our simulations are 
vertically resolved, with H/h at least 2. 

To ensure that the radial resolution is high enough in our sim- 
ulations, we resimulated the inner region of the disc employing 
the same number of particles used for our "fiducial" run. Figure 7 
presents a comparison between these high resolution runs, namely 
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Figure 8. Surface density normalized to initial value at two different times for simulations cold-adiabatic (red solid line), cold-iso (light blue dotted-dashed 
line) and for the analytical model (blue dotted line). Comparing the adiabatic with the isothermal case, both the position and the shape of the peaks are changed, 
as well as their height. 
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Figure 9. Surface density normalized to initial value at two different times for simulations hot-adiabatic (red solid line), hot-iso (light blue dotted-dashed line) 
and for the analytical model (blue dotted line). Comparing the adiabatic with the isothermal case, both the position and the shape of the peaks are changed, as 
well as their height. 



cold-iso-hres and hot-iso-hres, and the two simulations at standard 
resolution, namely cold-iso and hot-iso. We plot surface density 
as a function of radius at time t — 24d. After this time, the outer- 
most peak moves further than 1000i? s , thus reaching the outer edge 
of the disc. We do not find significant differences between the two 
runs, even if there is a factor 5 of difference in radial resolution. We 
can then conclude that the effects described in the previous section 
are not due to a lack of resolution. 



4.3 Adiabatic simulations 

In order to study the dependence of the evolution of the disc on the 
thermal physics, we also run simulations using an adiabatic equa- 
tion of state. In this case the gas can heat, and we expect thus a 
potential effect in the way the disturbances travel in the disc. 

In Figure 8 we plot the surface density, normalized to the ini- 
tial value, as a function of radius at time t = 48d and t — 97d for 
the simulations cold-adiabatic (red solid line), cold-iso (light blue 
dotted-dashed line), and for the analytical model (blue dotted line), 
while Figure 9 plots the same quantity for simulations hot-adiabatic 
(red solid line), hot-iso (light blue dotted-dashed line) and for the 
analytical model (blue dotted line). It can be noticed that the both 



position and the shape of the peaks are changed with respect to 
the isothermal case, and also their height. In particular, the simula- 
tion cold-adiabatic shows a shift in the position of the peaks that is 
someway in the middle between the simulations cold-iso and hot- 
iso. To show this, we plot in Figure 10 the results for simulations 
cold-iso (green dotted line), cold-adiabatic (red solid line) and hot- 
iso (light blue dotted-dashed line), zooming on the location of the 
peaks. It can be noticed that also the height of the peaks is smaller 
compared to the analytical model, and again intermediate between 
the simulations cold-iso and hot-iso. 

In order to account for these differences, we inspected the tem- 
perature structure of the cold-adiabatic simulation. In figure 1 1 we 
plot the temperature of the disc as a function of radius (circles) at 
time t = 97d, computed averaging the temperature of each SPH 
particle. For reference we also included the initial temperatures of 
the discs for the cold case (blue solid line) and for the hot one (green 
solid line). As can be seen, the disc heats up due to the shocks in- 
duced by the epicyclic motion. The total heating is not enough to 
reach the temperature of the hot case at the position of the out- 
ermost peak. Thus the A4 C parameter introduced in section 4.1 is 
intermediate between the values it assumes for the cold and the hot 
cases. This is consistent with the shape and position of the peaks 
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Figure 10. Surface density normalized to initial value at two different times 
for simulations cold-iso (green dotted line), cold-adiabatic (red solid line) 
and hot-iso (light blue dotted-dashed line). The simulation cold-adiabatic 
shows a shift and a height of the peaks that is someway in the middle be- 
tween the other two simulations. 




Figure 11. Temperature as a function of radius for simulations cold- 
adiabatic (circles, computed by averaging the temperatures of individual 
SPH particles), hot-iso (green solid line) and cold-iso (blue solid line) at 
time t = 97d. The temperature of the cold-adiabatic simulation at the posi- 
tion of the outermost peak is in the middle between the other two cases, thus 
accounting for the observed differences in the surface density evolution. 



observed, confirming that M t is a good parameter to characterise 
the behaviour of the disc after the mass-loss. In the same spirit we 
can interpret also the behaviour of the hot-adiabatic simulation, that 
shows height, shape and positions of the peaks with stronger differ- 
ences compared to the analytical model, since it can heat to even 
higher temperatures. 



5 DISCUSSION AND CONCLUSIONS 

In this paper, we have investigated the evolution of the dynamics 
of an accretion disc after the massloss of two merging black holes. 
Under the assumption of a disc composed of test particles, we were 
able to derive an analytical model for the surface density evolution 
of the disc following the mass loss. The model predicts the forma- 
tion of sharp density peaks in the disc. Once formed, these peaks 



travel outwards. However, we showed that they are not a type of 
hydrodynamical wave, but rather follow only from the kinematics 
of the disc. We derived also a timescale for the formation of these 
peaks, finding it of the same order of magnitude as previous find- 
ings. However, our model is fully time-dependent. 

To test the validity of our model, we set up numerical simula- 
tions, capable of taking into account the full hydrodynamics of the 
gaseous disc, using the SPH code PHANTOM. We found a good 
agreement in the shape and position of the peaks between the model 
and the simulations. There are however small differences, depend- 
ing on the disc parameters, in the position and shape of the peaks 
from the model predictions. To account for these discrepancies, we 
introduced the _M £ parameter. We showed that lower values of this 
parameter correspond to discs that deviate more strongly from the 
analytical model, because hydrodynamical effects are more impor- 
tant. We also showed that the timescale introduced is effective in 
predicting the formation of the density peaks. 

In the fluid disc, after some of the density peaks have formed, 
the epicyclic oscillations dissipate, thus causing the inner disc, 
where the dynamical timescales are faster, to significantly differ 
from the analytical model. We showed that the M t parameter is 
also effective in explaining how fast this dissipation proceeds, that 
is, how many peaks form in the disc. 

In our work we did not consider the radiation emitted by the 
fluid. We modelled two limiting cases, using an isothermal and 
an adiabatic equation of state, accounting for a very fast cooling 
timescale and a very long one, respectively. Further work is needed 
to quantify in which physical regime real discs lie. 
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